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I""! In the classical formulation of the Boussinesq approximation centrifugal buoyancy ef- 

^^ fects related to differential rotation, as well as strong vortices in the flow, are neglected. 

■ O However, these may play an important role in rapidly rotating flows, such as in astrophys- 

^ ical and geophysical applications, and also in turbulent convection. We here provide a 

'^-H straightforward approach resulting in a Boussinesq-type approximation that consistently 

^ accounts for centrifugal effects. We further compare our new approach to the classical 

• ^H one in fluid flows confined between two differentially heated and rotating cylinders. The 

^»>, results justify the need of using the proposed approximation in rapidly rotating flows. 



•^ 1. Introduction 

■^ In 1903 Boussinesq observed that: "The variations of density can be ignored except 

ly-. where they are multiplied by the acceleration of gravity in the equation of motion for the 



vertical component of the velocity vector" ( Boussinesq| 1903). This simple approximation 



f^ has had far-reaching impact on many areas of fluid dynamics; it allows us to approximate 

^' flows with small density variations as incompressible whilst retaining the leading order 

• • effects due to the density variations. Moreover, it is of great importance both analytically 

. ^H and numerically as it eliminates certain waves and acoustic modes which are more chal- 

rS lenging to treat. Many fluid dynamics problems have taken advantage of Boussinesq-type 

Cd approximations, rendering in most cases successful results in good agreement with ex- 

periments. However, some problems feature important physics neglected in the classical 
Boussinesq-type formulation. This is the case of fluids subjected to fast rotation. The 
traditional approach to account for centrifugal buoyancy is solely based on the angular 
velocity of the rotating reference frame, namely, the centrifugal term in the Navier-Stokes 
equations is proportional to fl^ and acting in the radial direction. This approximation is 
easy to apply in problems with a distinguished rotating frame of reference. Nevertheless, 
when such a frame of reference cannot be uniquely identified, the previous approximation 
neglects significant centrifugal effects due to differential rotation or strong internal vortic- 
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ity, which are specially important in fast rotating flows. The increasing interest on these 
flows because of their industrial (e.g. cyclonic dust collectors or vortex chambers) and 



scientific (astrophysical and atmospherical turbulence) applications (see Elperin et al. 



1998 1 motivates the development of a new approximation, which we here undertake. It 
is based on the classical Boussinesq approximation but it includes additional physical 
effects stemming from the advection term in the Navier-Stokes equations. This new for- 
mulation allows it to accurately cast rapidly rotating flows with mild variations of density 
into an incompressible formulation. In section §2, we describe a systematic way to achieve 
this, and we provide two different and easy to implement ways to account for centrifugal 
buoyancy effects in rotating problems. 

In order to illustrate the effect of different ways to account for centrifugal buoyancy in 
a Boussinesq-type approximation, we study numerically the linear stability of an axially 
periodic Taylor-Couette system with a negative radial gradient of temperature. Apart 
from its intrinsic interest, this setting has been widely used to model both atmospheric 



(e.g., see Randriamampianina et aL|2006 ) and astrophysical flows (e.g. see Petersen et al. 



[2007 ) where the fluid reaches high rotational speeds. 

Our simulations show that the classical Boussinesq approximation is valid in a wide 
range of Re numbers. However, for flows with high angular velocities important centrifu- 
gal effects arise. Here even the linear behaviour of the problem is significantly different 
for both approximations, justifying the application of the new approximation to account 
for centrifugal effects in rapidly rotating flows. The rest of the paper is organised as 
follows. In section §3 we give a detailed description of the system as well as the governing 
equations of the problem and its linearisation. A brief description of the basic flow is also 
provided. Section §4 introduces the Petrov-Galerkin method implemented to discretize 
the equations. In §5 the linear stability of the system considering both ways to introduce 
the centrifugal buoyancy is compared. Various cases of interest are analysed. In §5.1 we 
consider fluid rotating as a solid body, whereas in §5.2 shear is introduced in the system. 
We study quasi-Keplerian rotation in §5.2.1 and a system rotating close to solid body 
subjected to weak shear in §5.2.2. Discussion and concluding remarks are given in section 
§6. 

2. Boussines-type approximation for the centrifugal term 

In rotating thermal convection or stratified fluids the Navier-Stokes-Boussinesq equa- 
tions are usually formulated in the rotating reference frame, with angular velocity vector 
J7. The momentum equation in this non-inertial reference frame includes four inertial 
body force terms ( Batchelor||1967[ ), also called d'Alambert forces: 



pidt + u • V)u = -Vp + V • cr + «f - pV* 

(2.1) 
— pA — paxr — 2pr2xu — pfix (Jlxr). 

Here — pA is the translation force due to the acceleration A of the origin of the rotating 
reference frame, —pa x r is the azimuthal force (also called Euler force) due to the 
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angular acceleration a = dn/dt, -~2pu x fl is the Coriolis force and —pfl x (17 x r) is 



the centrifugal force (all of them per unit volume). In (2.1), p, p and u are the density, 
pressure and velocity field of the fluid, r is the position vector of the fluid parcel, and $ 
is the gravitational potential, so — pV<i> is the gravitational force. The term pi accounts 
for additional body forces that may act on the fluid. For a Newtonian fluid the stress 
tensor ex reads 

cr = -pI + /i(Vu + Vu^) + AV-uI, (2.2) 

where I is the identity tensor, p is the dynamic viscosity, and A is the second viscosity. 

2.1. The Boussinesq approximation in a rotating reference frame 

In most applications (geophysical flows, laboratory experiments), fl is constant and A = 
0, so we shall here focus on the gravitational, Coriolis and centrifugal terms. Note that the 
remaining inertial forces could be treated analogously. In the Boussinesq approximation 
all fluid properties are treated as constant, except for the density, whose variations are 
considered only in the "relevant" terms. Density variations are assumed to be small: 
P = Po + p'j with po constant and p' /po ^ 1; the p' term usually includes the temperature 
dependence, density variations due to fluid density stratiflcation, density variations in 
a binary fluid with miscible species of different densities, etc. With this assumption the 
continuity equation reduces to V • u = and the fluid can be treated as incompressible. 



As a direct consequence the shear stress term in the momentum equation (2.1 ) simplifles 
to the vector Laplacian, i.e. V • cr = pW^u. 

Identifying the relevant terms in the momentum equation is a more delicate issue. Any 



term in (2.1) with a factor p splits into two terms, one with a factor po and the other 
with a factor p'. If a po term is not a gradient, it is the leading-order term, and the 
associated p' term may be neglected. If the po term is a gradient, it can be absorbed into 
the pressure gradient and does not play any dynamical role, and therefore the associated 
p' term must be retained in order to account for the associated force at leading order. 
This is exactly what happens with the gravitational term: — poV$ = W{—po^), which is 
absorbed into the pressure gradient term and we must retain the —p'V^ term to account 
for gravitational buoyancy. The same treatment must be applied to the translation and 
centrifugal terms, yielding the gradient terms 

- poA - pon X (n X r) = V( -pol^ X rp - p^A ■ r), (2.3) 

as well as — p'A and — pTi x (J7 x r), which must be also retained. 



The Po part of the remaining terms in equation (2.1) (so far, we have considered the 
gravitational, centrifugal and translational forces) are not gradients, so they are retained 
as leading order terms and the corresponding p' terms are neglected, leading to the 
Boussinesq approximation equations in the rotating reference frame: 

poidt + u • V)u ^ -Vp* + pV^u + pi- p'V$ 

(2.4) 
— p'A — pqOl X r — 2pQ^ X u — pTi x (fi x r). 
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where 

p* =p + Po<i>--po\nxr\^+paA-r, (2.5) 

together with the incompressibihty condition V • u = 0. Of course, supplementary equa- 
tions are often needed; for example, if p' depends on the temperature, an evolution 
equation for the temperature must be included. 

2.2. Formulation in the inertial frame 

In many cases the fluid container is not rotating at a given angular speed, but different 
parts may rotate independently. For example Taylor-Couette flows with stratification 
and/or heating, cylindrical containers with the lids rotating at different angular velocities. 



etc. In these flows, there is not a natural or unique angular velocity fl to use in (2.4 1 and 



it may be more convenient to write the governing equations in the laboratory reference 



frame. In [2.2.1 we derive the momentum equation in the laboratory frame but for the 
sake of simplicity we assume that the fluid container rotates with angular speed ft. In 
§2.2.2| we show how the formulation is easily extended to account for the general case 
where a unique rotating reference frame cannot be identified. 

2.2.1. Formulation in the inertial frame: container rotating at angular velocity Jl 



The laboratory frame is an inertial reference frame, so the four inertial terms in (2.1 1 
are absent, and the balance momentum equation is 

pidt + V • V)v = -Vp + pV^v - pV* + p{, (2.6) 

where we have used v for the velocity field in the inertial reference frame, to distinguish 
it from the velocity in the rotating frame u. In order to implement the Boussinesq ap- 
proximation, we could naively repeat the previous analysis; since the only term which 
is a gradient is the gravitational force — /9oV$, we end up with an equation containing 
only the gravitational buoyancy, and the centrifugal buoyancy is absent. This appears 
reasonable, because the governing equations do not contain the rotation frequency ft of 
the container. However, Jl appears in the boundary conditions for the velocity, so it must 
be taken into account by a careful analysis of the nonlinear advection term. The easiest 
way to do this is by decomposing the velocity field as v = u -I- il x r, so the il x r part 
accounts for the boundary conditions (rotating container); u is precisely the velocity of 
the fluid in the rotating reference frame, with zero velocity boundary conditions. The 
advection term splits into four parts; 

V Vv = u • Vu + u • V{fl X r) + (J7 X r) • Vu + (J7 X r) • V{n x r). (2.7) 

Using the incompressibihty character of u, the dependence of fl on time but not on the 
spatial coordinates, and some vector identities, we can transform the advection term into 

V Vv == u • Vu + 2 rj X u + n X (O X r) + V X (u X (n X r)) . (2.8) 

We have recovered the Coriolis and centrifugal terms, and because J7 x ($7 x r) is a 
gradient, we must add a centrifugal contribution also in the inertial reference frame. 
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The last term in (2.8) accounts for the difference between the time derivatives in 



the inertial and rotating reference frames respectively. An easy way to see this is by 
considering the simple case where the two reference frames have the same origin, and ft = 
Jlk, where k is the vertical unit vector and fl is constant. Using cylindrical coordinates 
{r,9,z), with z in the vertical direction, we obtain 



V X (u X (O X r)) = fldeu. 
The change of coordinates between the inertial and rotating frame is 

r = r' , z = z', I 

e = e' + nt, t = t', J 



(2.9) 



(2.10) 



where {r',9',z') are the cylindrical coordinates in the rotating frame of the same fluid 
parcel with coordinates (r, 9, z) in the inertial frame; t and t' are the times in both 



reference frames. From (2.10) we obtain dt' = dt+VLdg, so the last term in (2.8), combined 



with dfVi results in the term SfU in the rotating frame. Finally, 9tv in the inertial 
frame contains an extra term, 9f(Jl x r) = a x r. Therefore, we have recovered all the 



inertial forces in the rotating frame (2.1 ), except for the translation force — pA, because 



in the example considered, (2.10), both reference frames have the same origin, and the 



translation is absent; by including a translation term in (2.10) we could also recover 



it. Now, the two formulations, including centrifugal buoyancy in both reference frames 
(inertial and rotating), fully agree. 

In the inertial reference frame, we are interested in a formulation in terms of the 
velocity field in the inertial frame v, instead of u as in (2.8). The analysis presented 



above considering the advection term results simply in an additional term, the centrifugal 
buoyancy. We have also discussed the effect of the decomposition v — u + Jlxrin 
the time derivative term. Now it only remains to consider the viscous term. However, 
V^(J1 X r) = because 17 x r is linear in the spatial coordinates and so its Laplacian is 
zero. The traditional Boussinesq approximation equations in the inertial reference frame 
are 



Po{dt 



V)v 



where p* = p + po'^ — 5Po|f^ x r 
V • u = 0. 



Vp* + pV^v + p f - /?' V$ - p'n X (fj X r), (2.11) 

^, and together with the incomprcssibility condition 



2.2.2. Formulation in the inertial frame: generalization 

We have shown that centrifugal buoyancy enters the governing equations via the bound- 
ary conditions and the advection term; no other term is affected in the Boussinesq ap- 
proximation. This now suggests a very simple formulation, consisting in keeping the 
whole density, p = po + p'., in the advection term. This formulation is easy to implement, 
and since most time-evolution codes for incompressible flows are semi-implicit (i.e. the 
viscous term is treated implicitly, whereas the advection term is treated explicitly), the 
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speed and efficiency of the codes do not change. The formulation reads 

pa{dt + V • V)v = -Vp* + AiV^v + pf - /9'V$ - p'(v • V)v, (2.12) 

where p* = p + po^, and ahows one to easily handle situations where different parts 
of a fluid container rotate independently. In these flows there is not a natural or unique 



angular velocity ft to use for a rotating reference frame in the formulation (2.11 ); however 



the angular velocities of the problem still enter the governing equations through the 



boundary conditions and the advection term. Hence formulation (2.12) provides a natural 



way to account for centrifugal buoyancy effects of these rotating flows in the inertial 
(laboratory) reference frame. This formulation is also appropriate if additional equations 
appear coupled with the Navier-Stokes equations, for example for large density variations 
in stratified flows. The treatment of the centrifugal effects can be carried out exactly in 
the same way presented here. 

2.2.3. Alternative formulation in the inertial frame and physical interpretation 



The extra term included in (2.12), p'(v • V)v, can be expressed in a different way, 



providing a closer resemblance to the expression in (2.11). Close to a rotating wall, the 
velocity field is v « Jl x r; this expression is exact at the wall (no slip boundary condition 
at a rigid rotating wall). The dominant part of the advection term is then 



(v • V)v w (J7 X r) • V{n X r) = J7 X (O X r) = -V(-|J7 x r|^) w -V( -v^). (2.13) 
As the dominant term is a gradient, it is necessary to include the p' term in the Boussinesq 



approximation. Replacing p'(v • V)v by —p''V{^v'^) gives the alternative form for (2.12) 
po{dt +v • V)v = -Wp 



1 



/iV^v + pf - p'V$ + p'V( -V- 



(2.14) 



This centrifugal effect is not only important when we have rotating walls, but also if a 
strong vortex appears dynamically in the interior of the domain; therefore it is advisable 
to always include this term in the Boussinesq approximation in order to account for all 
possible sources of centrifugal instability. 



We have presented two different ways, (2.12) and (2.14), of including the centrifugal 



buoyancy in rotating problems. One may wonder if there exists a canonical way to extract 
from the advection term the part that is a gradient, and then multiply this gradient by 
p'. The Helmholtz decomposition (Arfken & Weber 2005), writing a given vector field 



as the sum of a gradient and a curl, could serve this purpose, but unfortunately this 
decomposition is not unique (it depends on the boundary conditions satisfied by the curl 
part), and moreover it is not a local decomposition (i.e., in order to extract the gradient 
part, we need to solve a Laplace equation with Neumann boundary conditions). The two 



formulations presented here, (2.12) and (2.14), are simple and easy to implement, and 



deciding between one or the other is a matter of taste. 



The extra term we have included in (2.14), p'V(|v^), has an important physical in- 



terpretation; it is a source of vorticity due to density variations and centrifugal effects. 
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Taking the curl of (2.14) and using 



V X (v • Vv) = V X (a; X v) = V • Va; — a; • Vv, 
where a; = V x v is the vorticity field, results in an equation for the vorticity: 
Poidt + V • V)u = pqU} ■ Vv + ^V^a; + V x (pf) - Vp' x V$ + Vp' x V( -v^ 



(2.15) 



(2.16) 



The first three terms in the right-hand-side of (2.16) provide the classical vorticity evo- 



lution equation for an incompressible flow with constant density. The last two terms are 
the explicit generation of vorticity due to the gravitational and centrifugal buoyancies, 
respectively. In the following sections, we compare from the point of view of the linear 



stability the differences between the classical Boussinesq approach (2.111 and the first 



formulation proposed in this section, i.e. equation (2.12) 



3. Description of the system 

3.1. Governing equations 

We consider the motion of a fluid of kinematic viscosity v contained in the annular gap 
between two concentric cylinders of radii r^ and r,,. The cylinders rotate at independent 
angular speeds fli and flo- A negative radial gradient of temperature is also considered 
by setting the temperature of the inner cylinder to Ti = Tc + AT/2 and the outer cylinder 
to To = Tc — AT/2 where Tc is the mean temperature. We assume periodicity in the axial 
direction. The centrifugal buoyancy in the stationary frame of reference is included as in 



§2, equation(2.12): 

Poidt 



V • V)v = -Vp* -I- /iV^v - p'V$ - p'v • Vv, 



(3.1) 



where p* includes part of the gravitational potential, po*&- 

We asume p = pq + p' = po{l — aT), where T is the deviation of the temperature 
with respect to the mean temperature Tc, and po is the density of the fluid at Tc- We 
assume the gravity acceleration is vertical and uniform, so the gravitational potential 
is given by $ = gz; cylindrical coordinates {r,6,z) are used. With these assumptions, 
— p'V$ = poagTz where z is the unit vector in the axial direction z and a is the 
coefficient of volume expansion. The governing equations, including the temperature and 
incompressibility condition, are: 



(dt + v- V)v = -Vp + i^V^v + agTz + aTw ■ Vv, 
(9t -I- V • V)T = kV^T, 
V • v == 0, 



(3.2a) 
(3.2&) 
(3.2c) 



where k is the thermal diffusivity of the fluid. The equations are made dimensionless using 
the gap width d = Tq — r^ as the length scale, the viscous time <P /v as the time scale, 
AT as the temperature scale, and [y/d)^ for the pressure. In doing so, six independent 
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dimensionless numbers appear: 



Grashof number 
relative density variation 
Prandtl number 
radius ratio 

inner Reynolds number 
outer Reynolds number 



where Ap is the density variation associated with a temperature change of AT. In this 
system the Froude number is not particularly useful because we have two different rota- 
tion rates, fli and Qq, so the Froude number definition is not unique. 

From now on, only dimensionless variables and parameters will be used. The dimen- 
sionless governing equations are: 



G ^ agATd^i^^, 


(3.3a) 


e = a AT = Ap/po, 


(3.36) 


a — ly/n, 


(3.3c) 


V = n/ro, 


(3.3d) 


Rci = Qir,i,d/v, 


(3.3e) 


Rco = Vtorod/v. 


(3.3/) 



(a* -h V • v)v = 


-Vp + V^v + GTz + eTv • Vv, 


(3.4a) 


(9t + V • V)T = 


a-^V^T, 


(3.46) 


V • V = 0, 




(3.4c) 



The only change needed to recover the traditional Boussinesq approximation is to replace 
the centrifugal term eTv • Vv in (3.4a) by — efi^Trf, where f is the unit vector in the 
radial direction r. 

Due to the large number of parameters governing the system, a comprehensive explo- 
ration of the parameter space is extremely expensive, so we will vary Rei and G while 
keeping the remaining parameters constant (or a simple function of the two control pa- 
rameters Ri and G). In the present study, we keep fixed rj = 0.71, a typical value in 
experimental facilities, and a = 7.16, corresponding to water. 



3.2. Basic flow 

An analytical solution for the basic flow can be found by assuming only radial dependence 
for the variables of the problem. We also use the zero axial mass flux condition to fix the 
axial pressure gradient, i.e.: 



rwh(r)dr — 0. (3-5) 
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The resulting steady basic flow is given by: 



B 

r 



Ub{r) = 
Vb(r) = Ar 

Wb{r) = cl^Cir' - rf) + (^(r^ - rf ) + ^(r^ - r^)) 
„ / s 1 Inf r/r,) 



ln(r/rj 
Inr; 



p(r,z) 



Po 



ln?7 
^(4^+ 2 -4'- 



dr 



)z+ / (l--6Tfc(r))«2(^)_r 



(3.6a) 
(3.6&) 

(3.6c) 
(3.6rf) 
(3.6e) 



where Vb is the azimuthal velocity for the classical Taylor-Couette problem (Chan 



drasekhar 19611, whereas Wb and T{, correspond to convection in a conductive regime 



and appeared for the first time in Choi & Korpela (1980). The pressure varies linearly 



with the axial coordinate z, but the pressure gradient depends only on r, and therefore 
it is periodic in the axial direction. This axial pressure gradient mimics the presence of 



far away endwalls in any real situation, by unforcing the zero mass flux constraint (3.5) 



It is possible to give an explicit closed expression for p by integrating (3.6e), but it is 



quite involved and it does not appear in the problem solution. The expressions for the 
parameters A, B and C are: 

Rco - rjRci Rci - rjEco 



A 



C = - 



1 



B 



V: 



V '(i-'?)(i 

41nr7+(l-r;2)(3-772) 



r,')' 



16(1 - ?72)((l + 7^2) inry + 1 - 772) 



(3.7) 
(3.8) 



where (3.7 1 define the pure rotational flow in the azimuthal coordinate and C gives the 



axial component of the velocity field. The non-dimensional radii of the cylindrical walls 
are given by r^ = 77/(1 — 77), Tq = 1/(1 — 77). Note that the presence of the new centrifugal 
buoyancy term, proportional to e, does not modify the basic flow's velocity field, but only 
its pressure. 



3.3. Linearized equations 

We perturb the basic flow with infinitesimal perturbations which vary periodically in the 
axial and azimuthal directions, 



v(r, 9, z, t) = Vb(r) + e^("«+'=^)+^*u(r), 
r(r, 6*, z, t) = Tb{r) + e*(»«+'=^)+^*r'(7.), 



(3.9a) 
(3.9&) 



where Vb = {0,Vb,Wb) and Tb{r) correspond to the basic flow (3.6); u(r) = {ur,ug,Uz) 



and T'(r) are the velocity and temperature perturbations, respectively. The boundary 
conditions for both u and T' are homogeneous: u{ri) = u(ro) = T'{ri) = T'{ro) — 0. 
The axial wavenumber k and the azimuthal mode n define the shape of the disturbance. 
The parameter A is complex. Its real part A^ is the perturbation's growth rate, which 
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is zero at critical values, and its imaginary part Aj is the oscillation frequency of the 
perturbation. 



Using the decomposition (3.9) in the equations (3.4 1 and neglecting high-order terms, 



we obtain an eigenvalue problem, with eigenvalue A. It reads 



^ Id, dur. ru"^ + 1 

AUr = -^-[r^-j - Ur[ 5 

r or or r^ 



r 






2 

b ri^l 



^ 2 
2^1 



Id due ,n , 

r or or r^ 



r 



,dvb Vb 2in 

or r r^ 

Xu^^^-^ir^)-u,['^ + k^ + ^i'^ + kwb)il-en)] 
r or or r-^ r 



dwt 
dr 



(eTfe-lK+GT', 



or or or ct r^ r or 



(3.10a) 

(3.10&) 

(3.10c) 
(3.10d) 



Note that here the continuity equations and pressure terms are omitted because the 
Petrov-Galerkin method chosen to solve the resulting system of equations automatically 
satisfies the continuity equation and eliminates the pressure by using a proper projection 
(see next section). 



From (3.101 the equations for the traditional Boussinesq approximation can be easily 



obtained from (3.10) by setting e = in all terms except for —e{v^/r)T' . The traditional 
approximation contemplates only one rotating frame of reference for the system; the 



expression (3.6&) for the basic flow azimuthal velocity Wfc(r) — Ar + B/r has two terms, 
Ar corresponding to solid body rotation, and B/r corresponding to shear. It is natural 
to identify A as the frequency of the rotating frame of reference, fi^- In fact, if we take 
fli = il.o = rj, the Couette flow profile is: 



Vb{r) = Ar + 



B rtn 



n,rf 



(a, -r!o)(r,ro)2l 



fir = flrr, 



(3.11) 



and we recover the linearized version of the centrifugal term consisered by the classical 
approach, —efl^T'rr. In the general case with fli ^ Oq the traditional Boussinesq approx- 
imation is defined in the frame of reference rotating with fir = A. This approximation 
takes only into account the centrifugal buoyancy acting in the radial direction, which is 
obviously its main contribution. However, as we will see in Sj5l for high rotation rates 
other terms acting both in the radial and azimuthal directions become important and 
change the behaviour of the system. Part of the discrepancy stems from the fact that the 
effect of differential rotation is entirely neglected in the classical approximation. 
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4. Numerical method 

In order to solve numerically the eigenvalue problem described in the previous section, 
a spatial discretization of the domain must be made. This is accomplished by projecting 



the equations (3.10) onto a basis carefully chosen to simplify the process, 

^3 - {v e {^2{r^,ro)f I V • V = 0, v(r,) = v(ro) = 0}, 



(4.1) 



where {^2{i"i,ro))^ is the Hilbert space of square integrable vectorial functions defined 
on the interval [ri,ro), with the inner product 



v,u = 



V* • u rdr, 



(4.2) 



where * denotes the complex conjugate. For any v e V3 and any function p, we have 
(v, Vp) = 0. This consideration allows us to eliminate the pressure from the equations 
as we project them onto the basis. Moreover, the continuity equation is satisfied by 
definition of the space V3. For the temperature perturbation the appropriate space is 



Vi = {fe ^■2iu,ro) I fir,) = /(r„) = 0}. 
We expand the variables of the problem as follows 



X 



u(r) 
T'ir) 



^a,X, X, eVsxV^i, 



(4.3) 



(4.4) 



and projecting (3.10) onto V3 x Vi, we arrive at a linear system of equations for the 
coefficients Oj. 

The solution of the system is performed by means of a Petrov-Galerkin scheme, where 
the basis used in the expansion is different from the one used in the projection. The 
bases are composed of functions built on Chebyshev polynomials satisfying the boundary 
conditions. A detailed description of the method as well as the basis and functions used 



for the velocity field can be found in Meseguer & Marques (2000) and Meseguer et al. 



(2007), respectively. The basis functions for the temperature (last component of Xj in 



4.4), and for the projection (with ~) are: 



h,{r) = (1 - 2/2)T,_i(y), h,{r) = r\l - y^)T,^i{v), 



(4.5) 



where y = 2(r — r^) — 1 and Tj are the Chebyshev polynomials. As a result of this process, 
we obtain a generalized eigenvalue system of the form 



XMix = M2X, 



(4.6) 



where x is a vector containing the complex spectral coefficents(aj) and the matrices 
All and M2 depend on the parameters of the problem, the axial wavenumber k and 
the azimuthal mode n. This system is solved by using LAPACK. The numerical code 
written to perform this work implements the described method and analyses a range of 
fc, n and G provided by the user for a fixed Re number, searching for the critical values 
(3f?A = Ar = 0). Up to M = 200 radial modes have been used in order to ensure the 
spectral convergence when high Re numbers are considered. The code has been tested by 
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Figure 1 . Critical Grashof number Gc as function of inner cylinder Reynolds number Rei for 
fluid rotating as a solid-body. The solid line is the linear stability curve using the new approxi- 
mation for the centrifugal bouyancy proposed in this paper, whereas the dashed corresponds to 
the traditional Boussinesq approach. Different symbols indicate the two distinct mechanisms of 
instability. Up and down triangles represent the critical points due to the mechanism at moder- 
ate Rei for the new and traditional approximations respectively, whereas squares and diamonds 
correspond to the mechanism at large Rci. 



computing critical values for several cases in McFadden et al. ( 1984 ) and Ali & Weidman 



(1990), obtaining an excellent agreement with their results. 



5. Stability of differentially heated fluid between co-rotating cylinders 

We present a detailed comparison of the linear stability of the system using the tradi- 



tional Boussinesq approximation and our new approximation (3.10). We consider three 
different cases, all with 77 = 0.71 and a = 7.16. In the first one the cylinders are rotat- 
ing at same angular speed, corresponding to fluid rotating as a solid-body, and in the 
second and third cases the stability of a differentially rotating fluid, of importance in 
astrophysical flows, is considered in the presence of weak and strong shear. 

5.1. Cylinders rotating at same angular speed 
In this case a rotating frame of reference is readily identified and the shear term B/r in 



the basic flow azimuthal velocity ( 3.66[ ) is zero, whereas the term A corresponds to the 
angular velocity of the cylinders. Figure [l] shows the critical values of G as the rotation 
speed, indicated here by the inner cylinder Reynolds number Rci, is increased. In the 
case of stationary cylinders instability sets in at G = 8087.42, with kc — —2.74 and 
n — 0. The emerging pattern is characterized by pairs of counter-rotating torodial rolls, 
that unlike Taylor vortices have a non-zero phase velocity that causes a slow drift of 
the cellular pattern upward. Extensive information about natural convection instabilities 



can be found in the literature: Choi & Korpela (1980) and McFadden et al. (1984) for 



infinite geometries, and de Vahl Davis & Thomas ( 1969 ) and Lee et al. ( 1982 ) for finite 
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(a) (6) (c) 






Figure 2. Contours of the temperature disturbance T' at a z-constant section corresponding to 
the points marked as blue circles in figurejl] (a): Ret = 5- 10^ Gc = 21206.53. (6): Re, = 5.7- 10^, 
Gc = 33768.37. (c): Rd = 5.2 ■ 10^ Gc = 79670.16. There are 10 positive (dark gray; red in 
the online version) and 10 negative (light gray; yellow in the online version) linearly spaced 
contours. In all cases the critical azimuthal mode is n = 1 and k = O(10~'^). 

geometries. Obviously, without rotation both traditional (dashed line) and new (solid 
line) yield identical results. As rotation is increased no differences in the linear behaviour 
of the system are observed up to Rci ^ 5 x 10'"', where the two curves start to depart 
from each other. Up to this point and after a small initial region where several azimuthal 
modes up to n = 6 are involved, the basic flow loses stability to an azimuthal mode 
n = 1 with small axial wavenumber k ^ 10"'^. The shape of the critical modes along 
the stability curve is illustrated in figure [2J showing contours of constant temperature 
in a horizontal cross-section. The three states correspond to the circles in figure [l] and 
depict the transition between the lower and intermediate branches as we consider the new 
approximation. As we proceed forward along the critical curve the cold fluid progressively 
penetrates into the warm one and vice versa. The same behaviour is observed when the 
traditional approximation is used, nevertheless the values of Rci and Gc required are 
larger. 

As Rci increases beyond 5 x 10''^ the new terms in our approximation start becoming 
important and lead to different behaviour in the linear stability of the system. An analysis 
of the magnitude of each term in our approximation reveals that the differences observed 
in figure[l]at high Rci are due to terms involving the product VbUg, implying the existence 
of an important centrifugal force acting in azimuthal direction as high rotational speeds 
are reached. This provides evidence that the traditional formulation, including only the 
main (radial) contribution of centrifugal buoyancy, is a very good approximation if slow 
rotation is involved but other contributions may not be neglected in rapidly rotating 
fluids. Once the critical values given by both approximations differ, we can identify two 
interesting regions in parameter space. For Rci G [5-10^, 7.7-10^] the classical Boussinesq 
approach yields larger critical G than our approximation, whereas for Rci > 7 ■ 10^ the 
upper branch of the new approximation yields much lower critical values. Moreover, the 
differences keep increasing as Rci grows. 

The analysis performed reveals the existence of two mechanisms of instability associ- 
ated with the lower-intermediate and upper branches in figure [T] Different symbols are 
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Figure 3. (a) Critical axial wavenumber kc and (6) spiral angle of the modes arctan(fcc/w) as 
a function of Rei for the curves in figure [T] The inset is a close up at low Rei where the first 
mechanism stops being dominant and is superseded by spiral modes with angle far from 90 
degrees, indeed 0, corresponding to n = 0, at Rci = 0. 



used to represent the critical values corresponding to each mechanism in each problem. 
The differences between them are illustrated in figure |3J showing the evolution of the 
critical axial wavenumber kc and the angle of the spiral modes arctan(^^) versus Rci. 
Two regions with distinct characteristics are well-defined. The first mechanism of insta- 
bility has already been presented (see figure^. Low azimuthal wavenumbers, primarily 
n = 1, and very small axial wavenumbers characterize it. This corresponds to quasi two- 
dimensional modes and can be readily seen in figure IsFfc), showing that the angle of the 
spiral modes remains constant at about 90 degrees. The inset shows the small initial 
region where the spiral angle increases progressively until it reaches a vertical position. 
The second mechanism is characterized by n > 80 and kc ^ 0(1), also corresponding to 
quasi two-dimensional modes (see figure pp). Another common feature between the two 
types of instabilities is that the rotational frequency coincides with the angular veloc- 
ity of the container in both mechanisms and both approximations. This is in agreement 



with Maretzke et al. (2013), who have analytically proven that two-dimensional modes 



with k — Q always rotate at speed A (3.6&) in Taylor-Couette flows without heating. 



An interesting distinct feature of the second instability mechanism is localization near 
the inner cylinder. An example of these wall convection modes is shown in figure HI the 
critical disturbances are clearly different in the traditional and in our new Boussinesq 
approximations. 



5.2. Differentially rotating cylinders 
The traditional approximation for the centrifugal buoyancy neglects the part of the basic 



flow containing shear, i.e. the B/r term in (3.66). To quantify the influence of including 



shear in the centrifugal terms, we perform the same analysis as in the previous section 
but for differentially rotating cylinders. The amount of shear introduced is characterized 
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Figure 4. Contours of the critical disturbance temperature on a z-constant section for 
Rei — 700000. (a) Critical mode using the classical Boussinesq approximation. Here 
Gc = 489371.47, fee ~ 1-81 and n = 150. (b) Critical mode using our new approximation. 
Here Gc = 207906.92, fee = 0.54 and n — 116. In both cases, only 1/20 of the domain is shown. 
There are 10 positive (darker gray; red online) and negative (light gray; yellow online) contours. 

by the ratio of angular velocities /? = fti/flo; the further /3 is from unity, the stronger is 
the shear effect considered. 

5.2.1. Weak shear: rotation close to solid body (/3 — fli/flo = 1.006J 

We first consider the case where the container is rotating near solid body. Although 
shear may be here expected to play only a secondary role, this case serves the purpose 
of illustrating the importance of including shear effects in the centrifugal term. Figure [5] 
shows the neutral stability curve for the two approaches considered. Unlike the solid 
body case, the critical values Gc increase monotonically as Rei grows. This gives an idea 
of how important differential rotation is , even if the rotational speeds of each cylinder 
are only slightly different. Besides shear, centrifugal effects are also important in this 
configuration. From Rei > 2 x 10^ on the linear stability curves obtained by using both 
approximations become quite different. Similar features respect to the solid body case 
may be identified. At first the traditional approximation gives lower critical value of 
the Grashof. However, this region is smaller than in the solid-body case and ends at 
Rci ~ 2.9 X 10^ where both curves intersect. From that point on, the stability region 
predicted by the new approximation is smaller; the differences between the critical values 
given by both approximations keep increasing as larger Rei are considered. At the point 
where both curves first depart from each other Rei has half the value of that of the 
solid-body case. Consequently, the rotational speeds for which the new approximation is 
advisable are quite smaller in the presence of differential rotation. 

Critical axial and azimuthal wavenumbers exhibit similar behaviour to the solid-body 
case and so they are not shown here. Two mechanisms of instability are also found. 
The first one embraces the region 2 x 10^ < Rci and is characterized by fee ~ and 
1 < n < 6. Modes are similar to those obtained for the first mechanism in the solid body 
case. A subtle difference can be nevertheless pointed out. In the solid body situation the 
temperature disturbances fill the whole annulus, whereas differential rotation confines 
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Figure 5. Critical Grashof number Gc as function of inner cylinder Reynolds number Kei 
for rotation near solid-body (/3 = 1.006). Different symbols refer to two distinct instability 
mechanisms as in figure fT] 



(c) 






Figure 6. Contours of the critical disturbance temperature on a z-constant section, (a) n = 1, 
Rti — 178125, Gc ~ 92987.56. (6)-(c) Comparison of the traditional (fe) and new (c) approx- 
imation at Rei — 285000 (near the crossover point in figure [5| showing 1/20 of the annulus. 
(b) Gc = 154864.79, kc = 0.24 and n = 30. (c) Gc = 156547.54, kc = 0.39 and n = 75. There 
are ten positive (dark gray; red online) and negative (light gray; yellow online) linearly spaced 
contours. 



the perturbation towards the central part (see figure pa) . The second mechanism also 
presents the same features as in the solid body case, high azimuthal modes and kc G 
[0.5, 1.5], but differences in the flow appear that deserve to be highlighted. In the classical 
approximation the dominant wall modes are located at the inner cylinder, as it occurs 
in the solid-body case (figure^). In contrast, using the new approximation changes the 
location of the dominant wall modes to the outer cylinder (figure pp) . In view of these 
results we can say that considering shear effects in the centrifugal term of the Navier- 
Stokes equations may be extremely important: not only regarding the linear stability 
boundary but also the shape and location of the critical modes. 
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5.2.2. Strong shear: Quasi-keplerian rotation (j3 = Q.i/^o = 1.58J 

If l/rj > /3 > 1 the angular velocity decreases outward but the angular momentum 
increases. These flows, known as quasi-Keplerian flows, are used as models to investigate 
the dynamics and stability of astrophysical accretion disks. Here we choose a typical value 
j3 = 1.58 and as in the previous sections consider a negative temperature gradient in the 
radial direction, as expected in accretion disks. Figure [7] shows the neutral stability curve 
for the two approximations considered, as well as entirely neglecting centrifugal effects 
(e = 0). Here the three curves overlap, implying that shear is the completely dominant 
mechanism in this regime. As in the weak shear situation studied in the previous section, 
Gc is also characterized by a monotonous increase as Rei increases. Surprisingly, shear 
has a very strong stabilizing effect in this problem: without shear the critical Grashof 
number is ten times smaller at Re^ = 10^ than in the quasi-Keplerian case. 

Depending on the Reynolds number two mechanisms of instability are again found. The 
first mechanism exhibits a similar flow structure to that observed in the previous case. It 
also occurs at low Rei and is localized in the central part of the annulus due to the action 
of differential rotation. Figure p[a) shows the countours of the disturbance temperature 
in an horizontal plane. In contrast to what happens in the weak shear situation, these 
modes present a clear 3D structure with kc ^ —1. Small azimuthal wavenumbers are 
involved in this mechanism, ranging from n = 1 to n = 6. More remarkable differences 
are found when analysing the second mechanism. High azimuthal modes n ^ 50 arise as 
this mechanism becomes dominant, but unlike the solid-body and weak shear situations, 
the azimuthal wavenumber decreases as Rei increases. The same behaviour is observed 
in the axial wavenumber, so that the spiral angle quickly converges to 90 degrees as 
observed in the previous sections. Figure l8F6) shows that the instability is characterized 
by convection wall modes localized at the outer cylinder, as in the case of weak shear 
using our new approximation. Nevertheless, in quasi-Keplerian flows the dominant modes 
are always localized at the outer cylinder regardless of how centrifugal terms enter the 
equations. 

6. Summary and discussion 

We have identified some weaknesses in how the Boussinesq formulation is typically used 
to account for centrifugal buoyancy in the Navier-Stokes equations. In particular, this 
classical approximation neglects the effects associated to differential rotation or strong 
internal vorcitity. This has motivated us to develop a new consistent Boussinesq-type 
approximation correcting this problem. It consists in keeping the whole density in the 
advection term of the Navier-Stokes equations and thus it is very easy to implement in 
an existent solver. The new approximation allows it to accurately treat situations with 
differential rotation or when strong vortices appear in the interior of the domain, which 
may cause important centrifugal effects even in flows without global rotation. The latter 
may be especially relevant in simulations at high Rayleigh numbers (as e.g. in the quest 
for the 'ultimate regime', Ahlers et aL|2009|. Thus we argue that our formulation for the 
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Figure 7. Critical Grashof number Gc as function of inner cylinder Reynolds number Rei for 
quasi-Keplerian rotation (/3 — 1.58). The three curves differ only by about 1% and hence are 
only distinguishible in the inset (close up). 

(&) 





Figure 8. Contours of the critical disturbance temperature on a z-constant section, (a) 
Ret = 11681.03 with Gc = 4.1268 x lO", fee = -1.05 and n = 1. (6) Rei = 584051.72 with 
Gc ~ 1.8511 X 10^, fcc = 8.21 and n = 38. Ten positive (dark gray; red online) and negative 
contours (light gray; yellow online) are displayed. Only 1/20 of the domain is shown in (6). 



centrifugal terms should be always implemented whenever the Boussinesq approximation 
is used. 

The relevance of our new approximation has been illustrated with a linear stability 
analysis of a Taylor-Couette system subjected to a negative radial gradient of temper- 
ature. Three different cases have been studied. First, we have considered the container 
rotating as solid-body, i.e. without differential rotation. The critical values obtained for 
both traditional and new approximation agree up to Rei ^ 5.5 • 10^, where discrep- 
ances become significant. Beyond this point the conductive basic flow loses stability to 
quasi two-dimensional wall modes (aligned with the axis of rotation, as expected from 
the Taylor-Proudman Theorem) localized at the inner cylinder. Note that the large 
discrepancy in critical Grashof numbers observed at Rci € [5 • 10^, 10^] between both 
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approximations makes it possible to test them, against laboratory experiments. For ex- 



ample, in the experiments from Paoletti & Lathrop (20111, the experimental apparatus 



can easily reach Re = 10^, and the Grashof numbers about 5 x 10^ can be obtained with 
temperature differences about half a degree. 

We have also considered the case in which the cylinders rotate at different angular 
speeds, thus introducing shear. For weak differential rotation shear and centrifugal buoy- 
ancy effects compete and the critical values obtained with both approximations differ 
from each other at lower Rci ^ 2 x 10*^. Moreover, the new approximation gives rise to 
wall modes located on the outer cylinder, whereas the traditional approach yields wall 
modes on the inner cylinder as in the solid-body case. In strongly sheared quasi-Keplerian 
flows shear is so dominant that centrifugal terms may be entirely neglected in the linear 
stability analysis (discrepances in Gc are below 1% regardless of how centrifugal terms 
enter the equations, if at all). Here the critical modes are always localized at the outer 
cylinder. Surprisingly, shear has here a very strong stabilizing effect; at least when in- 
flnitesimal disturbances are concerned. Note that in more realistic models of accretion 
disks nonlinear instabilities have been found in similar regimes by |Klahr fc Bodenheimer 



(2003). Finally, it is worth noting that testing results with differential rotation in the 
laboratory is nearly impossible without including axial endwall effects. The large Rei 
involved will necessarily trigger instabilities due to the nearly discontinuous angular ve- 



locity profile at the junction between axial endwalls and cylinders (Avila|[2"012). 
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